!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine source_term (mx_new, my_new, mz_new, bx_new, by_new, bz_new, en_new, po_new,                                            &
      ro_old, mx_old, my_old, mz_old, bx_old, by_old, bz_old, po_old, gx, gy, gz,                                                  &
      bx_resis, by_resis, bz_resis, en_resis, ch2, cp2, dx, dy, dz, dt, nx, ny, nz)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: dx, dy, dz, dt, ch2, cp2

   real(8), dimension(nx, ny, nz), intent(inout) :: mx_new, my_new, mz_new, bx_new, by_new, bz_new, en_new, po_new
   real(8), dimension(nx, ny, nz), intent(in) :: ro_old, mx_old, my_old, mz_old, bx_old, by_old, bz_old, po_old, gx, gy, gz
   real(8), dimension(nx, ny, nz), intent(in) :: bx_resis, by_resis, bz_resis, en_resis

   integer(4) :: ip1, jp1, kp1, im1, jm1, km1
   real(8) :: div_b, s_mx, s_my, s_mz, s_bx, s_by, s_bz, s_en, dx2, dy2, dz2

!----------------------------------------------------------------------------------------------------------------------------------|
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   dz2 = dz * 2.0d0

   do k = 2, nz - 1
      kp1 = k + 1
      km1 = k - 1
      do j = 2, ny - 1
         jp1 = j + 1
         jm1 = j - 1
         do i = 2, nx - 1
            ip1 = i + 1
            im1 = i - 1

            div_b = ((bx_old(ip1, j, k) - bx_old(im1, j, k)) / dx2 +                                                               &
                     (by_old(i, jp1, k) - by_old(i, jm1, k)) / dy2 +                                                               &
                     (bz_old(i, j, kp1) - bz_old(i, j, km1)) / dz2) * div_clean_flag
            s_mx = -div_b * bx_old(i, j, k) + ro_old(i, j, k) * gx(i, j, k)
            s_my = -div_b * by_old(i, j, k) + ro_old(i, j, k) * gy(i, j, k)
            s_mz = -div_b * bz_old(i, j, k) + ro_old(i, j, k) * gz(i, j, k)
            s_bx = bx_resis(i, j, k)
            s_by = by_resis(i, j, k)
            s_bz = bz_resis(i, j, k)
            s_en = -(po_old(ip1, j, k) - po_old(im1, j, k)) / dx2 * bx_old(i, j, k) * div_clean_flag -                             &
                    (po_old(i, jp1, k) - po_old(i, jm1, k)) / dy2 * by_old(i, j, k) * div_clean_flag -                             &
                    (po_old(i, j, kp1) - po_old(i, j, km1)) / dz2 * bz_old(i, j, k) * div_clean_flag +                             &
                    mx_old(i, j, k) * gx(i, j, k) + my_old(i, j, k) * gy(i, j, k) + mz_old(i, j, k) * gz(i, j, k) +                &
                    en_resis(i, j, k)

            mx_new(i, j, k) = mx_new(i, j, k) + s_mx * dt
            my_new(i, j, k) = my_new(i, j, k) + s_my * dt
            mz_new(i, j, k) = mz_new(i, j, k) + s_mz * dt
            bx_new(i, j, k) = bx_new(i, j, k) + s_bx * dt
            by_new(i, j, k) = by_new(i, j, k) + s_by * dt
            bz_new(i, j, k) = bz_new(i, j, k) + s_bz * dt
            en_new(i, j, k) = en_new(i, j, k) + s_en * dt
            po_new(i, j, k) = po_new(i, j, k) * exp (-dt * ch2 / cp2 * div_clean_flag)
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine source_term
